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ABSTRACT 



Familial hemiplegic migraine (FHM) is a rare subtype of migraine with aura. A 
mutation causing FHM type 3 (FHM3) has been identified in SCN1A encoding the 
Navl.l Na + channel. This genetic defect affects the inactivation gate. While the Na + 
tail currents following voltage steps are consistent with both hyperexcitability and 
hypoexcitability, in this computational study, we investigate functional consequences 
beyond these isolated events. Our extended Hodgkin-Huxley framework establishes 
a connection between genotype and cellular phenotype, i.e., the pathophysiological 
dynamics that spans over multiple time scales and is relevant to migraine with aura. 
In particular, we investigate the dynamical repertoire from normal spiking (mil- 
liseconds) to spreading depression and anoxic depolarization (tens of seconds) and 
show that FHM3 mutations render gray matter tissue more vulnerable to spreading 
depression despite opposing effects associated with action potential generation. We 
conclude that the classification in terms of hypoexcitability vs. hyperexcitability is 
too simple a scheme. Our mathematical analysis provides further basic insight into 
also previously discussed criticisms against this scheme based on psychophysical and 
clinical data. 
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INTRODUCTION 

Familial hemiplegic migraine (FHM) is a rare monogenic, autosomal dominantly 
inherited syndrome with hemiparesis during the aura phase of migraine. Three distinct 
genetic mutations for FHM have been identified, in the CACNA1A calcium channel gene 
(FHM1), in the ATP1A2 Na,K-ATPase gene (FHM2), and in the SCN1A sodium channel 
gene (FHM3). It has been proposed that all three phenotypes reflect hyperexcitability in 
the form of increased susceptibility for spreading depression (SD) (van den Maagdenberg 
et al, 2007; Pietrobon, 2010). However, the functional connection between the molecular 
findings and a facilitated generation of SD is unclear. 

To determine the electrophysiological consequences of such a genetic defect, we 
integrate a mutation of FHM3 into three types of computational models of neuronal 
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dynamics. This allows us to bridge the gap between genotype and phenotype. A similar 
approach was used by Clancy & Rudy (1999). We use a standard Hodgkin-Huxley model 
for action potentials (AP) (Hodgkin & Huxley, 1952) and a model of SD (Htibel, Scholl 
& Dahlem, 2014) to evaluate the change in the threshold of generating SD by tolerating 
various brief intervals of transient ischemic attacks. Moreover, we use a model for anoxic 
depolarization (AD) (Zandt et al, 2011) that is derived from a seizure model (Cressman et 
al, 2009; Cressman et al, 201 1 ) as a test of the robustness of our results. 

The paper is organized as follows. In the Methods we introduce three computational 
models and our method to incorporate measured tail currents in FHM3 (Dichgans et 
al, 2005; Vanmolkot et al, 2007) into the Hodgkin-Huxley framework. In the Results 
we present simulations and analysis of the wild-type and mutant models. We end with 
the Discussion where we focus on three topics: (i) the appropriateness of the terms 
hypoexcitable vs. hyperexcitable, (ii) the seemingly paradoxically increased susceptibility 
to SD in the mutant model if one considers the firing rate, a measure that is usually used to 
quantify slow neural dynamics, and (iii) the inadequate concept of a threshold as a quantity 
measured by a single value. 



All three models are based on Hodgkin-Huxley type dynamics with different degree of 
complexity from the classical model to a second generation with time-dependent ion 
concentrations. 

Hodgkin-Huxley model 

The Hodgkin-Huxley (HH) model is one of the most widely used computational models 
in neuroscience. It is a conductance-based neuron model (Hodgkin & Huxley, 1952) and 
consists of four differential equations describing the membrane potential V and three 
gating variables m, n and h that determine the conductances of potassium and sodium 
channels. The change in membrane potential is proportional to the current that is flowing 
across the membrane with the proportionality constant given by the capacitance of the 
membrane C m . The individual currents are modeled as the conductance gi of the respective 
channel times the driving force, which is given by the difference between the membrane 
potential and the respective ion's reversal potential Ej, where i G {K, Na, leak}. Note that 
the conductance gj for voltage-gated channels, i.e.,;' e {K,Na}, is given by the maximal 
conductance gj times the respective gating variables as introduced below. The model takes 
into account a sodium current % a +, a potassium current a leak current Ji ea k that is 
carried by unspecified ions, and an applied current Z app . 

dV 1 

-77 - - 7T ( J Na+ + + heak ~ hpp ) , ( 1 ) 



METHODS 



W- — gN a m 3 h (V - £ Na ) , 

= gK» 4 (V - E K ) , 
keak=gl(V-Ekak)- 



(2) 
(3) 
(4) 
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Table 1 Model parameters for the Hodgkin-Huxley model. 


Name 


Value & unit 


Description 


Cm 


1 [tF/cm 2 


Membrane capacitance 


gNa 


120 m/cm 2 


Max. sodium conductance 


gK 


36 m/cm 2 


Max. potassium leak conductance 


a 


0.3 m/cm 2 


Leak conductance 


^Na 


50 mV 


Sodium reversal potential 


£k 


-77 mV 


Potassium reversal potential 


£ leak 


-54.402 mV 


Leak reversal potential 



In the HH model the potassium current is modeled as a delayed rectifier current with 
activation gate n while the sodium current is described by a transient current with an 
activation gate m and an inactivation gate h. All gating variables are voltage dependent and 
are given by the following equations: 



At X x 



with (5) 



Xoo = — -^-z- and (6) 

Ux + Px 
1 

= 

Olx + Px 



forx 6 {n,m,h}. (7) 



Xoq describes the steady-state of the gating variables and r x is the time constant. 
The rate equations for a x and ji x are voltage-dependent and given by 

0.1(V + 40) 



~ m l-exp(-(V + 40)/10)' (8) 
^ m = 4exp(-(y + 65)/18), (9) 
0.01(^ + 55) 



(10) 

l-exp(-(V + 55)/10) 
^ n = 0.125exp(-(y + 65)/80), (11) 
a/ ! = 0.07exp(-(y+65)/20), (12) 

B h = 1 . (13) 

l + exp(-0.1(V+35)) 

This model is capable of producing action potentials in response to depolarizations of the 
membrane caused by an appropriate externally applied current J app . All model parameters 
that were used in the simulations of the HH model can be found in Table 1 . It is interesting 
to remark that trying to study the effect of the mutation in a reduced two-dimensional 
model in the phase plane did not lead to promising results because the mutation quickly 
led to bistability, which is consistant with our results of a prolonged plateau of action 
potential and early depolarization block in the form of bistability. 
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Spreading depression model 

The classical HH model neglects the time-dependency of ion concentrations caused by 
spiking dynamics. Ions accumulate very slowly but also progressively due to the fluxes 
across the neuronal membrane. Therefore, changes in concentrations become significant 
either in the course of many rapid action potentials or under metabolic stress with 
insufficient ion pump activity, such as during transient ischemic attacks. Hence both the 
onset of spiking and also the response to reduced ion pump activity are of interest. These 
can be modeled by the spreading depression model described in more detail by Htibel, 
Schdll&Dahlem(2014). 

This model is also based on HH dynamics, but uses several changes and extensions. 
Instead of an unspecified leak current, a combined Na + -K + -leak current is used. The 
equations for sodium and potassium currents, including a pump current Ip that is 
introduced below, therefore change to 

W - (gL + A m ' h ) ■ (V - £ Na) + Hp , (14) 

'k+ = (^+^»V(^-£k)-2/ p . (15) 

Furthermore, the SD model uses dynamic ion concentrations to be able to model the 
breakdown of the ion gradients that is observed during SD. The intracellular potassium 
concentration K; and extracellular potassium concentration K e are modeled explicitly 
as dynamical variables, while the intra- and extracellular sodium concentrations (Na,- 
and Na e ) are computed from the potassium concentration due to the constraint of 
electroneutrality 

-r 1 = --I K +, (16) 
at a>i 

dK e y 

~T~ — +/diff(K e ) (17) 

at co e 

Na^Naf-IQ + Kf , (18) 
Na e = — (Na^ 0) - Na;) + Na<, 0) . (19) 



CO, 



The factor y converts currents to ion fluxes and depends on the membrane surface A m and 
Faraday's constant F: 

y = ^, (20) 

F 

a>i and co e are constants describing the intra- and extracellular volume, respectively, and the 
buffer flux J^tf is 

/diff = -Fdiff(Kbath - Kg) . (21) 

An overview of all constants and the values that were used in the simulations can be found 

in Table 2. 
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Table 2 Model parameters for the SD model. 



Name 


Value & unit 


Description 




1 |iF/cm 2 


Membrane capacitance 


J 

&Na 


0.01/5 m/cm 


Sodium leak conductance 




100 m/cm 2 


Max. gated sodium conductance 




0.05 m/cm 2 


Potassium leak conductance 


«£ 


40 m/cm 2 


Max. gated potassium conductance 






IjVjO aUQIUlll UUllCCllLld-LUJll 


"Ma 


1 70 mM /l 


T( C c/"\HinTYi /-~(~tTir~f>nt^i'citir\'n 
IvjO atJU.lLll.il LUllCCllLld-Llwll 


K, 


130.99 mM/1 


ECS potassium concentration 




4mM/l 


ICS potassium concentration 


^Na 


39.74 mV 


Sodium reversal potential 




-92.94 mV 


Potassium reversal potential 


ft); 


2,160 [im 3 


Volume of ICS 


We 


720 [im 3 


Volume of ECS 


F 


96,485 C/Mol 


Faraday's constant 


A m 


922 nm 2 


Membrane surface 


y 


9.556e 6 tlm2 c Mo1 


Conversion factor 


p 


5.25 ^A/cm 2 


Max. pump current 


4> 


3/ms 


Gating timescale parameter 


f diff 


3.75e— 5/ms 


Diffusion parameter 


K bath 


4 mM/1 


Potassium bath concentration 



If ion concentrations are time-dependent, they actually change drastically during 
neuronal activity To still maintain homeostasis an ion pump has to be included that 
pumps Na + ions out of and K + ions into the cell at a 3/2 ratio. The pump current thus 
depends on the extracellular potassium and the intracellular sodium concentration. The 
pump is modeled according to Barreto & Cressman (201 1 ) 

/^(Nai.Ke) = p(l + exp P 5 ~ ^ Xj (1 + exp(5.5 - K,))" 1 , (22) 

with p being the pump current strength. Note that the pump current also shows up in the 
equations for Na + - and K + - currents (Eqs. (14) and (15)). 

As a result of the dynamic ion concentrations also the reversal potentials become 
dynamic 

26.64 

£ ion = — — m([ion] e /[ion]j). (23) 

■Zion 

The fast gating dynamics of the m-gate is modeled adiabatically as 
m = moo (V). (24) 
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Note that in this model shifted versions of the rate equations are used (Cressman et al, 
2009; Cressman etal, 2011) 

0-1(7 + 30) 

a m = , (25) 

l-exp(-(V + 30)/10) 

y 0 m = 4exp(-(V + 55)/18), (26) 
0-01(7 + 34) 

a n = , (27) 

l-exp(-(V + 34)/10) 

y 6„ = 0.125exp(-(V + 44)/80), (28) 
a/ ! = 0.07exp(-(y+44)/20), (29) 

B h = 1 . (30) 

1+exp (-0.1(7 +14)) 

Furthermore, the time constants are scaled by a factor </> 

?x= 1 ■ (31) 

In contrast to Hubel, Scholl & Dahlem (2014) we did not reduce the dimension of the 
model further by assuming a linear or sigmoidal relation between n and h. Instead, h was 
kept dynamic since the changes caused by the mutation affect the /z-gate. 

Anoxia model 

As a test of the robustness of our results we investigate the effects of FHM3 also in a mutant 
model of anoxia (Zandt et al., 201 1 ) . In fact, migraine with aura has been linked to a higher 
risk of ischemic stroke (Kurth & Diener, 2012) . For furthere details on the rationale, see the 
Sec. Results. 

The anoxia model is similar to the SD model, but uses five more dynamic variables, in 
particular, it also models chloride ion dynamics. The other dimensions are due to explicitly 
modeling intra- and extracellular ion concentrations and not assuming mass conservation, 
and also electroneutrality is not assumed in this model. 

Therefore, in addition to Na + - and K + -currents as in Eqs. (14) and (15) a chloride 
(CI - ) channel is included, which contributes to the leak current 

dV 1 

-jt = -7^(W+'k++Jci) (32) 
dt L m 

Icr=S l a(y-E a )- (33) 

Intra- and extracellular ion concentrations are dynamic and modeled as 
dNa,- A 



W (34) 
dt VF Na 

dNa e BA , 
~ = — / Na + 35 
dt VF Na 

dCl; A , 

z = I n - 36 

dt VF 
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Table 3 Model parameters for anoxia model. 



Name 


Value & unit 


Description 




1 [J.F/cm 2 


Membrane capacitance 


*>JNa 


0.0175 m/cm 2 


Sodium leak conductance 


0 s 
«Na 


1 HO m /rin^ 

iuu in/ cm 


Jvlax. gated sodium conductance 


; 

ft 


0.05 m/cm 2 


Potassium leak conductance 


A 


40 m/cm 2 


Max. gated potassium conductance 




0.05 m/cm 2 


Chloride leak conductance 


Na; 


27 mM/1 


ECS sodium concentration 


Na e 


120 mM/1 


ICS sodium concentration 


K, 


130.99 mM/1 


ECS potassium concentration 




4 mM/1 


ICS potassium concentration 


-'-'IN a 


39.74 mV 


Sodium reversal potential 




-92.94 mV 


Potassium reversal potential 


<t> 


3/ms 


Gating timescale parameter 


A/VF 


0.044 SM/(S4) 


Conversion factor 




2.0 


Ratio ICS/ECS 


P 


28.1 |iA/cm 2 


Na-K-Pump rate 


G 


66 mM/s 


Glial buffering rate for K + 


6 


1.3 s" 1 


Diffusion rate 




4.0 mM 


Concentration K + in blood 


T 


310 K 


Absolute temperature 



dCl e 0A , 

~ - — J n - 37 

dt VF a 

dKi A , 

z = J K + 38 

dt VF ^ 

dK e pA , , 

dt VF 8 " 

The same pump current as in the SD model is used (Eq. (22)). While the total amount of 
sodium and chloride is constant, the extracellular potassium concentration can be buffered 
by glial cells (I g ) and diffuse into and out of the blood (J^) 

J g = G(l + exp(^^)) , (40) 

I d = e(K e -k OQ ), (41) 

h and n are dynamic and given by Eqs. (5), (6) and (25)— (31). The sodium activation gate 
m is adiabatically modeled as in Eq. (24) . For parameter values see Table 3 . 

Under physiological conditions this model behaves normally, as it responds with a single 
action potential to a short current pulse and with periodic firing when a larger current 
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-120 -80 -40 0 

V/ mV 



Figure 1 Inactivation time constant as a function of membrane potential. (A) Voltage-dependent time 
constant for mutation (r£) and wild- type (t^). Insets show the response of h to a voltage-clamp protocol. 
(B) shows the deinactivation (i.e., recovery from inactivation) process as a response to a step in voltage 
from —10 mV to —120 mV. (C) shows the inactivation process by stepping the voltage from —120 mV to 
— 10 mV. The intersections of the /j-curve with the 1 /e- and (1 — 1/ e)-lines, respectively, show the actual 
time constants. For deinactivation r£ is three-fold smaller than r^. For inactivation t ; * is three-fold lager 
than Tjj. 



of 1.5 mA/cm 2 or more is injected (not shown). This model is also able to show seizure 
activity (Cressman et al., 2009; Cressman etal., 2011). 

Modified time constant function based on tail currents 

The three models introduced above are given in their 'wild-type' formulation. The 
'mutant' formulation has only a single difference, a modified iNa current, as described 
in the following and illustrated in Fig. 1. 

From experimental data we know that the mutation leads to a two- to four-fold faster 
deinactivation (Dichgans et al, 2005) and to a two- to four-fold slower inactivation 
(Vanmolkot et al., 2007). We checked the robustness of our simulations within this range. 
The simulations presented here, however, were performed at an intermediate value of a 
three-fold change. 

To change the responsiveness of inactivation and deinactivation accordingly, we need 
to modify the time constant T/, of the gating variable h. In the mutant model this time 
constant is replaced by 

t£(V) = r h (V) ■ (/ci • tanh(cr • (V - V max )) + k 2 ). (42) 

The parameter V max shifts the sigmoidal tanh-function to the position of the maximum 
of the time constant function T/,(V). The slope factor of the sigmoidal tanh-function is 
er — 0.1 to ensure sufficiently rapid convergence to the limit of a three-fold change. The 
other parameters are K\ — 1.335 and K2 = 1.665. These parameters result from the two 
constraints k \ + K2 = f and K2 — K\ — 1 // for an /"-fold change. We chose/ — 3. 
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To test the mutant time constant r£, we simulated the experimental protocol performed 
by Dichgans et al. (2005) in the computational model. The membrane voltage is clamped 
to a holding potential of —120 mV and then stepped to a potential of —10 mV. At 
— 120 mV the /z-gate is completely deinactivated, i.e., open. The step to —10 mV causes 
the /z-gate to inactivate. Therefore, we can measure the time constant of inactivation with 
this protocol (see Fig. IB). In contrast, holding the membrane potential at —10 mV and 
then stepping back to — 120 mV allows us to measure the time constant of the process of 
deinactivation. At —10 mV the /z-gate is completely inactivated, i.e., closed, and the step 
to —120 mV causes the gate to deinactivate again, i.e., the gate reopens. An illustration 
of this protocol can be found in Fig. 1C. By using this procedure and measuring the 
two different time constants, it was assured that the chosen parameters lead to a 3 -fold 
slower inactivation and a 3-fold faster deinactivation. The main part of Fig. 1 shows the 
inactivation time constants T/, (black line) and tf (green line) for the wild-type and mutant 
model, respectively, as a function of the membrane potential V. 

Note that in the Hodgkin-Huxley formalism, the gating subunits of a channel are 
assumed to be identical and the inactivation and deinactivation as being independent. 
Therefore this formalism cannot represent certain dependencies in a straightforward 
manner in the kinetic states. For example, the inactivation of the Na + channel (represented 
by the /Vsubunit) has a greater probability of occurring when all subunits are open, 
therefore the inactivation depends on activation (represented by the three m-subunits). 
This violates the assumption of independent gating. Because of this independence in the 
HH formulation, the dynamics of the /z-gate is only described by a single time constant 
function T/,. An alternative ansatz is to use a Markov model to compute the occupancy of 
the channel in its various kinetic states as done by Clancy &Rudy (1999). 

RESULTS 

Three different models are investigated, a model of action potentials (AP), a model of 
spreading depression (SD), and a model of anoxic depolarization (AD). These models 
describe normal cell functions in terms of the dynamic repertoire either without genetic 
defect (three wild-type models) or with altered cell functions in FHM3 (three mutant 
models). The three mutant models (AP, SD, and AD) are the same as the wild-type models 
except that the Z>i a current has a different voltage-gating mechanism in the fast gating 
variable h. This is described in the wild-type model by the time constant T/, and in the 
mutant model by (see Methods). The observed functional consequences of FHM3 occur 
on time scales ranging from milliseconds to several tens of seconds. 

Mutant AP with marked plateau, increased responsiveness, 
delayed excitation block, and firing onset unchanged 

We first consider the shape of APs. The AP is rather directly affected by FHM3 through 
altered voltage gating in h. In other words, the results are consistent with the measured tail 
currents and therefore the results for a mutant AP are even to some degree predictable. This 
situation will change, when we model dynamics separated three orders of magnitude from 
AP dynamics. 
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I»PP / HA cm 2 

Figure 2 Spiking model. Comparison of wild-type (A) and mutant (D) spiking behavior. The main 
plots show bifurcation diagrams by varying the external current J app . For the wild-type model Hopf 
bifurcations can be found at J app = 9-78 |iA cm -2 and _T app = 154.52 |xA cm -2 . For the mutant model 
Hopf bifurcations occur at J app = 9.72 \iA cm -2 and J app = 175.02 |xA cm -2 . (B) and (E) show behavior 
in the oscillatory regime as a response to a constant input current of 12 \iA cm -2 . (C) and (F) show the 
response of the models in the excitatory regime to a 3 ms long current pulse with amplitude 3 |j,A cm -2 . 



For a single AP stimulated by a transient applied current / a pp(0 of 3 ms duration and 
3 (lAcm -2 amplitude (labeled 'excitatory' in Fig. 2), we observe that the mutant model 
compared to wild-type model leads to a prolonged AP with a marked plateau. This is 
consistent with the larger inactivation time constant ( V" e P) of the mutant as compared 
to the wild-type inactivation time scale r/,(V p ), cf. tail currents in Fig. IB. Note that we 
omitted before the explicit voltage dependency of the time constants, but now we make 
the dependency explicit because the mutant time constant function ^^(V) is in FHM3 
increased only for the regime of the membrane potential V being depolarised. This voltage 
regime is indicated by the superscript "dep" and it corresponds to an inactivation of h 
(closed h gates). 

Furthermore and a bit more subtle to observe, the mutant dynamics reacts faster 
to a sudden brief stimulation. The mutant model fires an AP that reaches its maximal 
amplitude just below 2 ms after the J app is turned off again, while in the wild-type model 
the maximal amplitude is reached only after about 3 ms. Again, this is also consistent with 
the defect in the time constant function r£ ( V). In this case it is explained by the decreased 
and therefore faster regime T/,(VT o1 ) compared to the wild-type. The mutant time constant 
function t£(V) is decreased for V being in the polarised resting state indicated by the 
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200 




Figure 3 Nonlinear firing-rate function F(I a pp) for wild-type model (black, solid) and mutant model 
(green, dashed). 



superscript "pol", cf. tail currents in Fig. 1A. This is the regime of deinactivation (open h 
gates). 

The modified AP profile is also observed during spiking, i.e., in the oscillatory regime, 
when a constant / app larger than — by definition (see below) — the rheobase current / r h 
is applied. Individual APs in the spike train show this plateau (labeled 'oscillatory' in 
Fig. 2). As a result the spiking frequency is reduced in the mutant model, despite the overall 
increased responsiveness (Fig. 3). This decreased spiking frequency can be associated with 
hypoexcitability as the neural response is usually characterized by the firing-rate function. 

To get some further quantitative measures of the effects of FHM3 with regard to 
excitability, we investigated the change of stability in the resting state by varying the 
input current J app . This is a bifurcation analysis (Fig. 2). The determined two so-called 
bifurcation points mark the beginning and end of the oscillatory spiking regime. The 
first Hopf bifurcation point (HB1) is the onset of oscillation at a minimal value of / app , 
which is the definition of the rheobase current J r h- For the wild-type model the first Hopf 
bifurcation (HB1) is at/™ 1 = J rh — 9.78 \iAcm~ 2 and the second Hopf bifurcation (HB2) 
at J™ 2 = 154.5 ^iA cm -2 , which determines the excitation block as the oscillation ceases at 
this point. For the mutant model these Hopf bifurcations occur at J r h — 9.72 \iA cm -2 and 
I™ 2 = 175.0 ^lAcm -2 . The first Hopf bifurcations (HB1) are subcritical, while HB2 are 
both supercritical. This means that if the / app is not slowly ramped towards the rheobase 
current J^, one can observe the oscillatory regime even before the 1^. Hence the two 
firing-rate functions in Fig. 3 start slightly before the values given here for HB1, with the 
mutant model starting again earlier. 

With regard to the rheobase current, the values for the wild-type and mutant differ 
by less then 0.6%, with the mutant value being smaller, which, at least in principal, 
corresponds to hyperexcitability, though due to the small magnitude this seems negligible 
for all practical purposes. However, the excitation block observed at the second critical 
transition HB2 occurs at larger values of / app for the mutant model. The mutant channels 
tolerate an increased maximal /™ 2 by 13% compared to the wild-type. This means that 
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the mutant neurons exhibit oscillatory behavior in a larger range of applied currents. 
Therefore, this shift establishes a gain-of- function, which indicates hyperexcitability. 

To summarize, while the reduced firing frequency indicates hypoexcitability, increased 
responsiveness and delayed excitation block indicate hyperexcitability 

Mutant more vulnerable to SD 

We now focus on effects of FHM3 upon cellular functioning that occurs in the same neural 
substrate that generates APs but on time scales at least three orders of magnitude separated 
from AP dynamics, that is, effects that occur during several tens of seconds up to minutes. 
This is the time scale of SD. It is therefore relevant for pathological conditions, for instance, 
in migraine with aura. In accordance with this pathophysiological context, we select the 
stimulations of SD in the wild-type and mutant model as rather large perturbations to 
neural homeostasis such as a compromised energy supply during focal hypoperfusion 
that induces and occurs in conjunction with migraine aura symptoms (Olesen et al., 1993; 
Friberg et al, 1994). 

In particular, we investigate the effect of a breakdown of the Na + -K + -pump upon the 
membrane potential V and reversal potentials E^ a and £k- For this purpose the maximum 
pump rate p is linearly down-regulated to 20% of its physiological value within 10 s, then 
p is kept at 20% for a variable time window, and finally p is linearly up-regulated back to 
100% within 5 s. The stimulation trace of p is shown in Fig. 4 with the dashed-dotted line. 
The specific choice of the variable time window is additionally marked for the wild-type 
and mutant stimulation trace by an annotated two-headed arrow. Let us remark that in our 
studies we also used two other perturbations, namely a transient increase in extracellular 
K + concentration (by increasing Kbath) and a large current pulse / a pp, with basically the 
same results (not shown). 

We determined the minimal duration of the variable time window with reduced pump 
rate (20%) that is just no longer tolerated and results in a long lasting but transient 
breakdown of the reversal potentials i?Na and £k characteristic for SD. For this purpose 
we increased the variable time window by 0. 1 s steps. While the wild-type model could not 
tolerate a period of 13.6 s of reduced pump rate at 20%, the mutant model was less robust 
and could not tolerate a period of 7.2 s of reduced pump activity (Fig. 4). Therefore, the 
mutant model is approximately only half (53%) as robust to periods of reduced ion pump 
activity as the wild-type model is. 

Shorter stimulation periods did not lead to full blown SD signals. In this case, the 
spiking ceased about a second after the interval began that increased the pump rate 
back from 20% to 100% (this interval lasts 5 s) and, more importantly, both membrane 
potential V and reversal potentials E^ and £k recovered within only a few seconds back 
to physiological values (not shown). Thus, SD profiles of these potentials, which followed 
longer stimulation periods, are clearly distinguished by a all-or-none phenomenon. Not 
only do membrane potential V and reversal potentials E^a and £k change dramatically 
after the stimulation is off, but also full recovery from SD to the initial physiological values 
takes very long. Of course recovery reaches the resting state only asymptotically. For up to 
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Figure 4 Spreading depression model. Development of SD in wild-type (A) and mutant (B) models. A 
SD is elicited by down-regulating the pump current to 20% of its maximal value for 1 3.6 s (wild-type) and 
7.2 s (mutant), respectively (see blacked dashed line). The red and blue dashed lines show the temporal 
development of the sodium and potassium reversal potentials. 



one to two hours the changes in particular in E^a are observable, while the signals in Fig. 4 
are shown only for 100 s. It is noteworthy that the neuronal state is already back to basic 
functioning emitting APs if stimulated after the repolarization, that is, even if the resting 
state is not fully recovered. Similar dynamics is described in other computational models of 

SD by Kager, Wadman & Somjen (2000), Yao, Huang & Miura (201 1 ) and Hubel, Scholl & 
Dahlem (2014). 

To summarize, in terms of susceptibility to SD the mutant model is hyperexcitable. This 
seems to be in contrast to the major effect of the mutant upon the AP firing frequency that 
indicates that the mutant model is hypoexcitable (Fig. 3). This will be further discussed in 
the Discussion. 

Effects in anoxia model consistent with SD model 

Last, we study a model of AD (Zandt et ah, 2011); the AD model shares many features 
with the SD model but is more detailed (see Methods) and hence effects obtained with 
this model serve as control to compare them with effects obtained from the SD model. 
The model was first published to study slow waves after decapitation in a computational 
model (Zandt et ah, 2011). By repeating this with a mutant version of this model, our focus 
is set very similar to the previous section. In the decapitation study, anoxia is modeled by 
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Figure 5 Anoxia model. Response of wild-type (A) and mutant (B) membrane potential to a complete 
breakdown of pump, glial and diffusion currents at t = 5 s (black dashed vertical line). Red and blue 
dotted lines show the Na + and K + reversal potentials over time. The time from the onset of spiking until 
the beginning of the excitation block is approximately 6.7 s without and 2.7 s with mutation. 



completely switching off all pump, glial, and diffusion currents, see Fig. 5. In fact, Fig. 5A 
with the wild-type model is a reproduction of the simulations performed by Zandt et al. 
(2011). 

Note that patients with migraine with aura are at greater risk for stroke (Kurth & 
Diener, 2012). Thus there is a rationale to perform this comparison beyond the mere 
confirmation of plausibility of our results obtained above with the SD model. However, 
the multiplicity of potential links include not only common genetic risk factors but also 
indirect links like common triggers outside the brain, e.g., microemboli caused by cardiac 
shunts. Furthermore, the model investigated by Zandt et al. (201 1 ) is derived from a model 
suggested by Cressman et al. (2009). This model exhibits periodic bursting similar to 
seizure activity. Both migraine and epilepsy have genetically based forms caused by various 
mutations in genes, while the mutation in FHM3 differs markedly within the several 
mutations in SCN1A therein that it is not associated with epilepsy (see Introduction). 
Investigating the underlying ion homeostasis in the three conditions of epilepsy, migraine, 
and stroke may yield interesting results in future investigations of computational models 
that can unify certain dynamical aspects and link disease genotype to phenotype. However, 
this is clearly beyond the scope of this study. 
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In this study, let us only refer to the dynamics resulting from switching off all pump, 
glial, and diffusion currents until the excitation block and compare the wild-type and 
mutant model. After a gradual rise of the membrane potential that lasts in either case about 
30 s (note that the simulated 'decapitation' occurs in Fig. 5 at t = 5 s), the membrane 
potential reaches the AP threshold, subsequently resulting in a final burst of spiking. These 
initial, less than a minute lasting, phases in the wild-type and mutant model are indeed 
very similar to the initial phases in the SD model following a transient energy failure. A 
minor difference is that the gradual rise is overall slower, but this is explained by a slightly 
different geometry (larger extracellular space) and by the chloride ion dynamics (Hubel, 
Scholl & Dahlem, 2014). The similarity supports the robustness of our results, as this model 
is an established model showing anoxia (Zandt et al, 201 1 ) and seizure activity (Cressman 
etal, 2009; Cressman et al, 2011). 

To summarize, also for AD the slow gradual fall of the potentials does not significantly 
differ during the initial leak phase in the wild-type and mutant model, while once the 
model is spiking the excitation block occurs about 2.5-times faster, corresponding to a 
faster breakdown of ion gradients due to spiking, in the mutant model. 

DISCUSSION 

Our main result is that the mutant model is more susceptible to spreading depression (SD). 
With our computational model, we bridge the gap between the tail currents measured 
by Dichgans et al. (2005 ) and altered cell function that constitutes the phenotype of 
migraine with aura. Importantly, in a computational model we can follow in all needed 
detail how the complex interactions of channel dynamics lead to altered cell function. A 
similar approach was taken, for instance, to link a genetic defect to its cellular phenotype in 
a cardiac arrhythmia by Clancy &Rudy (1999). 

In the discussion, we mainly highlight aspects of hypoexcitable vs. hyperexcitable and 
the concept of a threshold. 

Hypoexcitable vs. hyperexcitable 

The increased susceptiblility to SD does not contradict the reduced firing frequency for a 
given stimulation current J app , although this change in firing frequency indicates that the 
mutant model is hypoexcitable. 

Firing a single action potential (AP) is a form of cellular excitability manifested as a 
transmembrane voltage jump without significant changes in ion concentrations. SD is a 
form of cellular excitability manifested by massive changes in ion concentrations. There is 
not necessarily a direct relation between the two excitable systems, not even with regard 
to merely classifying terms such as hypoexcitable and hyperexcitable. Rather, AP and SD 
can be viewed as largely independent phenomena, because while sharing the same neural 
substrate, AP and SD are separated by time scales differing in three orders of magnitude 
(see below). Notwithstanding, the massive breakdown of ion gradients in SD is, of course, 
mediated by APs that occur on the fast time scale. 

In our view, "hypoexcitable" vs. "hyperexcitable" is in any case too simple a classification 
scheme even considering AP and SD in isolation on their respective time scale. To support 
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this criticism of classifying neural dynamics in migraine, let us mention that this problem 
was also addressed in the psychophysical and clinical contexts, see studies by Shepherd 
(2001 ) and Coppola, Pierelli & Schoenen (2007) and references therein; further support 
comes from the mathematical picture (below) — which are two sides of the same coin. 

To illustrate this with only a single example, consider, as already mentioned above, 
that the mutant channels exhibit an increased range of spiking activity with a delayed 
excitation block by 13% compared to the wild-type. We argued that this larger spiking 
range establishes a gain-of- function. Consider further the increased responsiveness of the 
mutant model. Both indicate a form of hyperexcitability with regard to AP. In contrast, 
the change in firing frequency of AP indicates at the same time that the mutant model is 
hypoexcitable (Fig. 3). 

SD susceptibility 

How do these three diverse effects observed for APs (delayed excitation block, increased 
responsiveness, and lower firing frequency) manifest on the longer time scale under the 
condition of SD? 

In terms of susceptibility to SD, the shifted excitation block (see HB2 in Fig. 2) might 
misleadingly suggest that the mutant model is less susceptible to SD. This is similar to 
the lower firing frequency that we considered above. Since the characteristic sustained 
breakdown of the reversal potentials and £k is ignited in our model only if the system is 
driven by any stimulation into the excitation block, its delay in the mutant model seems to 
suggest that a longer stimulation may be needed and therefore a higher threshold exists. 

To show the actual situation in Fig. 4, we highlighted a critical time window by a gray 
shade. This critical time window opens with start of the reduced pump rate recovery (from 
20% back to 100%) and it closes with the beginning of the excitation block. Considering 
only the delay of the excitation block and the low frequency, it may seem surprising at 
first, this critical period lasts 3.4 s in the wild-type model and only 2.5 s in the mutant 
model. Note that this 'paradox' can also be observed in the overall shorter duration of 
the whole initial firing pattern in the mutant SD model. Our attention should be on 
signals that can actually be measured in a clinical setting, hence our focus is on these 
signals also in the presentation of the computational model, where we can "measure" 
everything. The reduced pump rate corresponds to hypoperfusion signals. The excitation 
block in SD corresponds to the first peak in an electroencephalography (EEG) signal, 
cf. the work by Zandt et al. (201 1 ) where the simulated membrane potential is averaged and 
high-pass filtered, cut-off at 0.1 Hz, to estimate the EEG — although this EEG might only be 
observable intracranially. 

That the mutant model is more susceptible to spreading depression (SD) is exclusively 
explained by the much larger amount of ions transferred across the membrane during 
spiking. This, in particular the intracellular ion concentration, cannot easily be measured 
even in an in vitro setup. In Fig. 4, we see this by the much steeper slope of the reversal 
potential £k in the mutant model. 
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The multidimensional concept of thresholds 

The complex question about susceptibility requires a deeper understanding of what a 
threshold is. In fact, the very reason why we have to get beyond the idea of "hypoexcitable" 
vs. "hyperexcitable" as a useful characterization of the system (see above) is that there is no 
one-dimensional ansatz to determine a threshold as a demarcation. 

Before explaining this further, let us give one more explicit example. In other model 
variants of SD (Kager, Wadman & Somjen, 2000), a stimulation of SD may even stop before 
the excitation block is reached. In this case a sustained afterdischarge carries the system 
into the depolarisation block that then marks the start of the actual SD events. Clearly, in 
this case the depolarisation block cannot be considered being the actual threshold, because 
the system is 'before' this point when the stimulation is already off again. 

In general, excitability or all-or-none phenomena do not possess a threshold in terms 
of single quantity, whether it is a particular membrane depolarisation that demarcates 
the all-or-none response in the form of an AP or a critical duration of hypoperfusion 
that demarcates the all-or-none response in form of SD. A detailed analysis of neural 
models shows that a threshold is a multidimensional surface (manifold) not a single 
number as first shown by FitzHugh (1955) and as discussed in a modern style by Mitry 
etal. (2013) and applied to migraine by Dahlem (2013). So the actual use of computational 
models goes far beyond numerical simulations. We gain a deeper understanding of the 
principal mechanisms in precise mathematical relationships, of which we can only give a 
very general overview in this paper. 
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